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SUMMARY 

A numerical  procedure  for  computing  the  laminar  flow-field  in  nozzles  at 
throat  Reynolds  numbers  of  300-3000  is  described.  Such  nozzles  are  found  in 
spacecraft  position  control  thrusters,  chemical  lasers,  and  low  density  hyper- 
sonic wind-tunnels.  A parabolic  approximation  of  the  Navier-Stokes  equations 
(the  'boundary  layer  equations')  is  transformed  to  Von  Mises  form  and  solved 
numerically  by  a central  difference  method.  The  entire  subsonic  and  supersonic 
flow-field  is  computed.  No  assumptions  or  approximations,  other  than  those 
inherent  in  the  flow  equations,  are  involved.  The  method  is  of  the  direct  type 
and  the  flow-field  for  any  given  nozzle  geometry  may,  in  principle,  be  computed. 
Examples  are  given  comparing  computed  results  with  published  experimental  data. 
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INTRODUCTION 


Viscous  nozzle  flows  occur  in  at  least  three  fields  of  technology.  In  low 

thrust  satellite  position  control  thrusters  the  nozzle  throat  size  is  of  order 

one  millimetre,  and  typically  the  throat  Reynolds  number  is  in  the  range 

1000-3000.  Usually  they  operate  with  a propellant  of  nitrogen,  hydrogen  or 

ammonia  or  mixtures  of  these  gases,  often  at  high  temperature.  In  low  density 

hypersonic  wind-tunnels  the  nozzles  are  larger,  having  throat  radii  from  several 

millimetres  to  several  tens  of  millimetres  and  operate  with  throat  Reynolds 

numbers  in  a similar  range,  usually  with  nitrogen  or  air  . The  nozzle  walls  are 

often  cooled.  In  chemical  lasers  the  nozzles  are  of  the  same  size  range  as 

satellite  control  thrusters  and  operate  in  the  same  Reynolds  number  range,  often 

2 

with  fluorine  as  the  'working  fluid'  . 

The  numerical  procedure  described  in  this  Report  was  developed  for 

predicting  the  performance  of  satellite  control  thrusters.  It  is,  however, 

equally  applicable  to  any  nozzle  flow,  although  in  this  Report  only  laminar 

flows  are  considered.  The  basic  drawback  to  operating  a convergent-divergent 

nozzle  at  laminar  Reynolds  numbers  is  that  viscous  shear  forces  in  the  flow  in 

the  divergent  part  of  the  nozzle  inhibit  the  acceleration  of  the  flow.  The 

boundary  layer  thus  thickens  and  in  extreme  cases  may  encompass  the  majority  of 

the  mass  flow.  In  the  case  of  propulsion  nozzles  this  results  in  poor  nozzle 

efficiency  and  a reduction  in  specific  impulse  (thrust  per  unit  mass  flow).  In 

the  case  of  low  density  hypersonic  wind-tunnel  nozzles  it  limits  the  obtainable 

Mach  number  and  limits  the  region  of  constant  Mach  number  in  the  flow  core.  In 

the  case  of  chemical  lasers  the  nozzle  flow  viscous  forces  reduce  the  stagnation 

pressure  available  for  recovery  in  the  laser  cavity,  and  this  has  a direct 

bearing  on  the  output  power.  In  addition  the  reactions  which  produce  the  laser 

power  occur  between  the  oxidiser  and  fuel  streams.  This  is  usually  limited  to 

2 

the  nozzle  boundary  layer  so  the  influence  of  the  boundary  layer  on  the 
operation  of  the  laser  is  important. 

Most  viscous  nozzle  flow-field  analyses  belong  to  one  of  two  types.  In 

one  type  of  analysis  the  flow  is  approximated  by  a viscous  boundary  layer  and  an 

inviscid  core.  The  flow  in  the  boundary  layer  may  be  obtained  by  using,  for 

3 

example,  Cohen  and  Reshotko's  approximate  integral  method  and  then  using  an 

iterative  technique  to  couple  the  boundary  layer  and  inviscid  core.  This 

4 5 

approach  was  followed  by  Murch  et  at  and  by  Potter  and  Carden  . Alternatively, 
the  equations  of  the  boundary  layer  may  be  transformed  by  introducing  a 
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stream-function  and  applying  a modification  of  the  Lees-Dorodnitsyn  transformation 
(a  combination  of  the  Mangier  transformation  and  Illingworth  transformation) . 

The  resulting  coupled  pair  of  partial  differential  equations  may  then  be 
solved  directly  by  a finite  difference  procedure  (Whitfield  ) or  the  equations 
may  be  simplified  and  transformed  to  a pair  of  ordinary  differential  equations 
by  assuming  local  similarity.  These  equations  may  then  be  solved  more  readily 
at  the  expense  of  the  constraints  imposed  by  similarity.  This  approach  was 
adopted  by  Edwards  . In  both  cases  the  boundary  layer  solutions  are  coupled 
iteratively  to  the  one-dimensional  core  flow. 

In  the  other  type  of  analysis  the  entire  flow-field  is  computed  by  solving 
the  Navier-Stokes  equations  directly.  Only  one  analysis  of  this  type  has  been 

g 

reported  . Here,  Rae  solved  a parabolic  approximation  of  the  Navier-Stokes 

9 

equations  (the  'boundary  layer  equations')  in  the  channel  flow  approximation 
for  the  whole  flow-field.  An  implicit  Crank-Nicholson  finite  difference  scheme 
was  used. 

At  low  Reynolds  numbers  experimental  evidence***  indicates  that  the 
inviscid  core  almost  disappears.  The  viscous  boundary  layer-invise  id  core 
approximation  is  then  not  even  approximately  true.  It  then  becomes  very 
desirable  to  account  for  the  viscous  effects  over  the  whole  flow,  and  an  analysis 
of  the  complete  flow-field  becomes  necessary.  It  is  assumed,  of  course,  that 
the  Reynolds  numbers  are  not  so  low  that  the  Navier-Stokes  equations  cease  to  be 
valid. 

A diff iculty  wi th  the  channel  flow  approximation  is  that  the  channel  walls  are 
limited  to  small  angles  of  convergence  and  divergence.  Furthermore,  no  account 
is  taken  of  any  radial  pressure  gradient.  In  general  the  streamtubes  are 
inclined  to  the  symmetry  axis.  There  will  thus  be  a component  of  the  pressure 
gradient  along  the  streamtubes  in  the  radial  direction,  and  experimental  data  in 
Ref  10  shows  that  this  component  may  not  be  ignored. 

The  numerical  analysis  reported  here  is  of  the  'entire  flow-field 
computation'  type.  The  nozzle  flow-field  is  computed  using  a form  of  the 
boundary  layer  equations  in  which  a radial  pressure  gradient  exists. 

2 THE  FLOW  EQUATIONS 

The  starting  point  is  the  form  of  the  boundary  layer  equations  developed 
by  Probstein  and  Elliott**.  These  equations  were  also  used  by  Whitfield*1  and 
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by  Edward*7.  Symbols  are  defined  in  the  nomenclature.  The  coordinate  system 
is  shown  in  Fig  I. 
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where  p is  the  density,  h is  the  static  enthalpy,  p the  static  pressure, 
u the  laminar  viscosity,  k the  thermal  conductivity,  C the  specific  heat 
at  constant  pressure,  u the  velocity  in  the  longitudinal  direction,  v the 
velocity  in  the  transverse  direction,  x the  longitudinal  coordinate,  y the 
transverse  coordinate,  and  r is  the  radius. 

1’robstein  and  Elliott  derived  these  equations  for  external  flow  where  the 
boundary  layer  thickness  is  much  smaller  than  the  longitudinal  radius  of 
curvature.  These  equations  are  here  extended  to  internal  flow  and  to  the  whole 
of  the  flow-field.  The  distance  x is  measured  from  some  reference  plane  along 
a given  streamtube  and  the  distance  y is  measured  from  the  symmetry  axis  along 
a surface  orthogonal  to  the  streamtubes.  In  this  case  this  is  also  a surface 
of  constant  static  pressure  since  the  transverse  momentum  equation  is  neglected 
leading  to  3p/3y  ■ 0 . In  this  approximation  this  arises  from  neglecting  the 
centrepetal  force  component  which  accompanies  a change  in  direction  of  the  flow, 
to  39/3x  is  assumed  negligible.  This  is  a good  approximation  in  the  divergent 
part  of  the  nozzle  where  the  streamtubes  are  only  slightly  curved,  but  in  the 
region  of  the  nozzle  throat  this  approximation  is  not  so  satisfactory. 

The  equations  are  thus  not  limited  to  small  angles  of  convergence  and 
divergence  of  the  nozzle  walls,  but  are  limited  to  small  first  derivatives.  In 

g 

this  sense  this  approximation  is  one  order  higher  than  that  used  by  Rae  and 
there  is  a radial  component  of  the  pressure  gradient. 
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2. 1 Transformation  of  the  flow  equations 

Equation  (3)  may  be  written  in  parametric  form: 


1 2 

where  the  parameter  «|i  is  the  compressible  stream  function  . From  the  chain 
rule  for  differentiation  we  have  the  following  operators: 


(to,  - to.  • to,  to. 


(to.  ■ to.  (to. 


so  that 


(&\  ‘ (&),  - pvr(^), 
‘ our(&)„  • 


Thus  equations  (I)  and  (2)  may  be  transformed  from  coordinates  x,  y to 
x,  <|>  and  we  have: 


3u  I dp  3 ( 2 3u\ 

37  - "^d7+37Vwpur  w) 


3h  1 dp  __  3 / 2 k 3h\  2/3uV 

37  * p d5 + 37  (,pur  c;jfj  + wpur  \fa)  • 


This  is  known  as  Von  Mises  transformation  . 

The  value  of  ^ at  the  axis  of  the  nozzle  is  taken  to  be  zero  (an 
arbitrary,  but  convenient,  choice).  At  the  nozzle  wall: 


Xw 

■/ 


purdy 
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and  is  independent  of  x since  the  total  mass  flow  is  conserved  throughout  the 
nozzle.  Hence: 


where  m is  the  total  mass  flow  rate. 

It  is  known  that  in  the  Von  Mises  form  the  boundary  layer  equations  are 

. ' . 14 

likely  to  lead  to  a more  stable  numerical  procedure  . In  addition  they  may  be 

solved  numerically  by  a marching  integration  procedure  which  avoids  any  iteration 

as  will  be  shown  later. 

This  method  of  handling  the  boundary  layer  equations  is  essentially  the 

same  as  that  employed  by  Patankar  and  Spalding  for  general  second  order  parabolic 

...  15 

equations 

Equations  (II)  and  (12)  are  written  in  terms  of  non-dimensional  variables 
(defined  in  the  ncmenc lature)  and  become: 

- ! ! — + 8Pr  fu*p*u*r*2  (15) 

3x*  u2  p*u*  dx*  rr0  V P 5*V  V J 


Y0  ' 1 


(p*u*r*2  K- 

p*  dx*  3iJ>*  Cp*  3i|>* J 


2 


where  M is  the  Mach  number,  y the  ratio  of  specific  heats  and  Pr  is  the 
Prandtl  number.  The  superscript  (*)  refers  to  non-dimensional  variables  and  the 
subscript  (0)  refers  to  mean  values  at  the  duct  entrance. 

2.2  The  boundary  conditions 

An  idea  of  the  nozzle  geometry  may  be  gained  from  Fig  3 or  Fig  6.  It  will 
be  seen  that  at  the  nozzle  inlet  the  wall  is  parallel  with  the  symmetry  axis. 

It  is  assumed  that  at  the  nozzle  inlet  the  flow  has  constant  enthalpy  across  its 
cross-section  and  that  a fully  established  (parabolic)  velocity  profile  exists. 
It  is  further  assumed  that  there  is  no  radial  velocity  component.  Thus  we  have: 


2(1  - **)■ 


1 at 
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At  the  nozzle  axis  the  conditions  express  the  symmetry.  So  we  have: 

3u^  AKA 

3^"°»  v*«0,  ~_.0  at  1//*  - 0 . (18) 

It  is  assumed  that  at  the  nozzle  wail  non-slip  conditions  are  appropriate, 
ie  velocity  slip  and  temperature  jump  are  ignored.  Hence: 


u*  ■ 0 , 


v* 


0 , 


3h*  cp*q* 

jp  “ at  **  - > 


(19) 


where  q*  is  the  non-dimensional  heat  flux  at  the  nozzle  wall. 

For  the  two  test  cases  presented  in  section  5 the  wall  was  assumed  to  be 
adiabatic, ie  q*  ■ 0 . Note  that  the  transverse  enthalpy  gradient  is  left  in 
terms  of  the  transverse  coordinate  y*  . The  reason  for  this  is  that  as 

y - yw  . 3h/3*  - - because  u - 0 and  it  is  possible  to  incorporate  the 

boundary  conditions  as  shown  and  thus  avoid  the  singularity. 

3 THE  GAS  PROPERTIES 

The  computations  presented  in  this  Report  were  performed  with  nitrogen. 
The  equation  of  state  used  is: 


p - pRT/m  (20) 

where  R - 8.3143  x |03  j k"1  kg-mole'1 
and  m * 28.013  (the  molecular  weight). 

The  gas  is  thus  assumed  to  be  ’ideal’. 

The  Prandtl  number  was  assumed  constant  at  Pr  - 0.7. 

The  viscosity  was  assumed  to  obey  Sutherland's  law: 

14  5 T3^  -7  o 

v - TH76  * t * ,0  N • (2i) 

The  specific  heat  is  assumed  constant  at  C ■ 1.0397  x 103  J kg”1  k”1 
and  the  ratio  of  specific  heats,  y - 1.400  . 

It  is  assumed  that  all  energy  states  are  in  mutual  equilibrium  and  frozen 
flow  losses  are  ignored. 
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4 THE  FINITE  DIFFERENCE  PROCEDURE 

Equations  (15)  and  (16),  together  with  the  initial  and  boundary  conditions 
(17)  to  (19),  the  equation  of  state  (20)  and  the  viscosity  law  (21)  and  various 
mainly  geometrical  relationships  (which  will  be  developed  in  section  4.3)  form 
the  starting  point  for  the  numerical  procedure. 

A central  difference  technique  is  used  to  solve  equations  (15)  and  (16). 

A general  point  in  the  mesh  scheme  and  the  domain  of  integration  in  the  x*  - if/* 
plane  is  shown  in  Fig  2.  The  equations  have  been  integrated  up  to  and  including 
the  mth  axial  point  and  the  solution  is  about  to  be  extended  to  the  m + 1th 
axial  point;  $ is  an  arbitrary  function  and  B is  an  arbitrary  coefficient. 

Thus  we  have: 


3» 
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, n)  l^m+ 1 , 


n+1 
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2&4i* 
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2A** 
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.14,16 


The  Crank-Nicolson  scheme  is  sometimes  not  stable  for  flow  in  a duct 
so  a fully  implicit  scheme  is  used  and  values  of  4>  are  evaluated  along  the 
line  m + 1.  Values  of  the  arbitrary  coefficient  B are  taken  along  the  m line 
where  the  previous  solution  is  known.  This  avoids  non-linearity  and  the 
resulting  necessity  for  iteration. 

Equations  (15)  and  (16)  may  be  written  in  the  form: 


9$ 

9x* 


A Si  + 3^  (B  w)  + c 


(25) 


for  4 = u*  we  have: 
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*oMo 


p*u* 


(26) 
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B = ®Pr0  W*P*u*r* 


C = 0 


(27) 

(28) 


for  $ = h* 
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Y0  ~ 1 1 
Y0  P* 


B = 8p*u*r* 


2 
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c = 


2,  .x  y*p*u*r* 


8Pr0VY0  " 0 


4 (Aip*) 


r (5- 


n+l 


u*  y 

m,  n-  IJ 


(29) 


(30) 


(31) 


Note  chat  in  equation  (16)  3u*/3<l>*  is  evaluated  along  the  m line: 


/ 3u*\2  1 / . * \2 

W 4(A^*)2  \m»n+1  m‘n'V 


(32) 


In  each  case  the  coefficients  A*  B and  C are  evaluated  on  the  basis  of 
the  known  solution  at  the  previous  step.  They  are  therefore  regarded  as 
constants. 

4. 1 The  flow  equations  in  finite  difference  form 

Writing  (25)  in  finite  difference  form  using  (22)  and  (24)  we  have: 


*m+l 


♦ ♦. 


o+l 


♦ ♦. 
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Then  we  have: 


23  24 
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\ \ \ 

\c  ^ c X c 

N-l ,N- I N-l ,N  N-l ,N+1 


fra+ 1 , 1 
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where 


A .Ap*  + Ax*C  . ♦ $ 
m,  l m.i  m.i 


Note  that  the  flow  has  been  divided  into  N equal  movements  of  stream  function: 


(37)  gives  a system  of  N-l  equations  in  N+l  unknowns. 

4 . 2 The  boundary  conditions 

The  boundary  conditions  give  zero  gradient  of  $ at  the  symmetry  axis. 

So  o and  c,,  may  be  added  to  c,,  , the  first  column  deleted 

m+ 1 , 1 m+ 1,2  11  J 12 

from  the  coefficient  matrix  in  (37),  and  $ . . deleted.  The  boundary  condi- 

m+ 1 , i 

tions  at  the  nozzle  wall  are  of  two  types.  Either  the  value  of  $ is  specified 

or  its  first  derivative.  When  $ = u*  the  first  of  these  alternatives  applies 

and  we  have  $ ■ 0 . In  this  case  the  last  column  may  be  deleted  from 

m+ 1 , 1 

the  coefficient  matrix  and  p , ..  . deleted.  When  p = h*  , and  an  adiabatic 

ID+  I , N ♦ 1 

wall  is  assumed  we  have  ♦ , so  that  c„  , „ . may  be  added  to 

m+l,N  m+l,N+l  N-l, N+l 

c^_j  jj  » c^e  last  column  deleted  from  the  equations  and  $m+|  deleted.  The 

effect  is  the  same  in  both  instances;  the  coefficient  matrix  is  reduced  to 
tridiagonal  form  and  we  have  N-l  equations  in  N-l  unknowns: 


t 


12 


(C11  + C.2> 
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: 23 


22  '23  24 
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\ \ \ 
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rm+l,2 


m+ 1 ,N 


N- 1 


(40) 


Other  conditions  at  the  nozzle  wall  may  be  treated  in  a similar  way.  If 

$ is  specified  and  is  other  than  zero  (for  example  the  wall  temperature  may  be 

given,  or  the  velocity,  if  a velocity  jump  is  allowed),  then  the  terra 

c»,  , may  be  subtracted  from  b„  . and  the  last  column  deleted  and 

T)-l,N+rm+l,N+l  1 N-l 

ip  deleted. 

m+1 ,N+1 

The  case  where  a finite  gradient  is  specified  is  a little  more  involved. 
For  example  suppose  the  heat  flux  is  given.  Then  at  the  nozzle  wall: 
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Hence 
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(44) 


In  finite  difference  form: 
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The  RHS  is  known  on  the  basis  of  the  previous  step.  So  h* 


raav 


be  expressed  in  terms  of  h*  . M and  the  coefficient  matrix  reduced  as  before. 

ID+  I ) N 


4.3  Solution  of  the  finite  difference  equations 


The  two  tridiagonal  systems  of  equations  for  u*  and  h*  are  solved  at 
each  step  by  a Gaussian  elimination  method' For  a tridiagonal  system  of 


equations  this  results  in  a very  compact  procedure. 

The  temperature  vector  is  computed  from  the  enthalpy  vector,  and  the 
density  vector  from  the  temperature  and  pressure. 

The  radial  coordinate  is  computed  from  a recurrence  relation: 


pur 


3y 


(46) 


so 


fJjT  " 2p*u*r* 


(47) 


It  is  easily  shown  from  the  geometry  that 


3r* 

3y* 


cos  9 


(48) 


so 


3i|)» 

3r* 


2p*u*r* 


cos  0 


(49) 


In  finite  difference  form: 


4A>1>*  cos  0 


Ar* 


r*  . - r* 
n+ 1 n 


n+l 


(p*  + p*  (u*  + u*  /r*  + r*  ,\ 
V n 'n+l;  \ n n+l)  v n n+l ) 


(50) 


Hence: 


r* 

n+l 


4A\|>*  cos  0 


n+l 


+ P 


n+l 


IK 


* 


(51) 
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Thus  given  that  r*  ■ 0 and  values  for  the  density,  velocity  and  stream- 
tube  inclination  at  the  previous  step  the  radii  may  be  easily  computed. 

In  general  the  final  radius  r*+J  does  not  quite  coincide  with  the  duct 

wall,  ie  the  computed  flow-field  does  not  quite  fit  the  duct.  This  is  because 

the  procedure  uses  the  pressure  and  pressure  gradient  computed  at  the  previous 

step.  The  procedure  corrects  the  pressure  to  make  the  computed  flow-field  fit 

the  duct  and  the  resulting  pressure  and  pressure  gradient  are  used  at  the  next 

axial  step.  Thus  no  iteration  is  required,  the  procedure  simply  marches  along 

the  nozzle,  the  pressure  and  pressure  gradient  being  controlled  by  an  error 

parameter.  The  necessary  correction  to  the  pressure  is  estimated  from  the  one- 

1 8 

dimensional  flow  equation  connecting  pressure  change  and  area  change  : 

. (52) 

P l - M2  A 

The  mean  Mach  number  across  the  flow  is  used  and  the  area  change  is 
computed  from  the  radius  error.  It  was  found  desirable  to  use  only  a fraction 
of  the  pressure  correction  to  ensure  stability.  In  addition  the  magnitude  of 
the  correction  must  be  limited  to  avoid  problems  with  the  singularity  at  M = 1. 

Because  the  radii  fluctuate  the  streamtube  inclination  also  fluctuates, 
and  this  is  coupled  back  into  the  equations  and  produces  instability.  In  order 
to  overcome  this  effect  a small  step  length  is  used  (1/250  of  the  throat  radius) 
and  the  change  in  inclination  is  computed  every  7-25  steps  depending  on  the  wall 
curvature.  In  regions  of  high  curvature  the  smaller  interval  was  used  since 
the  radius  is  changing  more  rapidly  and  small  errors  have  less  effect  on  the 
inclination.  The  effect  is  to  smooth  the  streamlines,  and  it  successfully 
stabilises  the  procedure  at  the  expense  of  greatly  reducing  its  efficiency. 

The  streamtube  inclination  0 is  calculated  from  the  geometry.  We  have: 

sin  0 « ~ (53) 

dX 

where  Ar  is  the  change  in  radius  from  the  chosen  point  7-25  steps  upstream, 
to  the  current  point  along  a given  streamtube.  Ax  is  the  corresponding  change 
in  longitudinal  distance.  The  expression  becomes  increasingly  inaccurate  as 
30/3x  increases.  It  can,  however,  be  shown  to  be  consistent  with  the  other 
assumptions,  specifically  3p/3iJ>  ■ 0 . 
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If  the  flow  parameters  are  such  that  the  flow  does  not  choke  at  the 
nozzle  throat  the  pressure  decreases  up  to  the  nozzle  throat,  and  then  increases 
the  flow  being  everywhere  subsonic.  If  the  procedure  tries  to  force  too  great 
a mass  flow  through  the  nozzle  it  fails;  the  flow  chokes  before  the  throat  is 
reached.  When  this  situation  occurs  the  error  between  the  computed  flow-field 
and  the  duct  radius  increases  rapidly,  resulting  in  a greatly  increased  negative 
pressure  gradient.  The  pressure  then  becomes  negative  and  the  procedure  fails. 
Both  these  situations  may  be  easily  detected  and  upper  and  lower  bounds  on  the 
mass  flow  rate  established.  The  procedure  then  takes  the  mean  of  the  upper  and 
lower  bounds  and  establishes  a new  upper  or  lower  bound  depending  whether  the 
new  mass  flow  rate  leads  to  premature  choking  or  no  choking.  Thus  the  flow  rate 
is  adjusted  until  the  flow  chokes  at  the  nozzle  throat  and  then  the  pressure 
decreases  monotonically  throughout  the  duct.  This  bisection  method  of  iteration 
has  only  linear  convergence,  but  convergence  is  guaranteed  and  usually  only 
about  eight  iterations  required.  This  point  is  explained  further  in  the  next 
paragraph. 

The  eigenstate  corresponding  to  choking  is  not  well  defined.  If  the  mass 

flow  rate  is  slightly  greater  than  the  'true'  eigenvalue  the  computed  flow-field 

expands  slightly  in  the  region  of  the  throat,  to  accommodate  it.  The  sonic 

surface  then  intersects  the  nozzle  axis  slightly  upstream  of  the  geometric  throat. 

If  the  mass  flow  is  very  carefully  adjusted  to  make  the  error  parameter  a 

minimum  in  the  throat  region  the  sonic  surface  intersects  the  axis  slightly 

downstream  of  the  geometric  throat.  It  is  thus  possible  to  make  the  sonic 

surface  intersect  the  axis  at  any  arbitrary  point  (within  limits)  around  the 

nozzle  throat.  It  is  possible  to  vary  the  mass  flow  by  about  0.5Z  at  a throat 

Reynolds  number  of  2000  and  achieve  a smooth  expansion  to  supersonic  flow.  At 

a Reynolds  number  of  500  the  flow  is  more  critical  and  it  is  possible  to  vary 

the  flow  only  by  about  0.05Z.  This  phenomenon  is  a result  of  the  numerical 

0 

method  used,  at  least  in  part.  Rae  also  noted  a similar  effect,  except  that 
his  mass  flow  was  much  better  defined  and  the  limits  between  which  the  sonic 
surface  could  be  varied  were  both  downstream  of  the  nozzle  throat.  The  effect 
is  very  convenient  from  the  practical  standpoint.  It  means  that  when  the  whole 
procedure  is  iterated  to  find  the  choking  eigenstate  it  is  only  necessary  to  find 
the  mass  flow  rate  to  three  significant  figures,  which  saves  time. 

Note  that  no  additional  assumptions  or  fitting  (as  was  done  in  Ref  8)  is 
required  in  the  region  of  the  singularity  at  the  throat.  When  the  choking  mass 
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flow  rate  has  been  determined  the  procedure  will  start  from  its  initial  condi- 
tions and  continue  to  compute  through  the  subsonic  flow  through  the  throat  and 
into  the  supersonic  flow  in  the  expansion  cone. 

The  computations  presented  in  this  Report  were  performed  in  single 
precision  - seven  significant  figures  - on  a large  minicomputer.  Double 
precision  made  no  significant  difference  to  the  results.  The  procedure  is  very 
inefficient  and  for  the  two  test  cases  presented  the  total  time  for  each  was 
about  4 hours.  This  is  a total  of  eight  iterations  of  the  flow  up  to  the  throat 
to  determine  the  eigenstate  and  one  complete  determinator  of  the  flow  to  the 
nozzle  exit.  This  time  is  reduced  to  a total  of  about  20  minutes  using  a more 
modern  computer  (PRIME  400).  Thus  although  procedures  such  as  the  one  described 
in  this  Report  are  very  lengthy  the  advent  of  modern  fast  computers  means  that 
they  become  more  practicable. 

5 RESULTS 

A survey  of  the  literature  reveals  little  detailed  experimental  data  on 
the  flow-fields  in  low  density  nozzle  flows.  The  exception  is  Ref  10.  Here, 

Rothe  measured  the  gas  density  and  temperature  at  a large  number  of  points 
throughout  the  flow  of  nitrogen  through  small  nozzles  using  an  electron  beam 
fluorescence  technique.  His  measurements  covered  throat  Reynolds  numbers  (based 
on  throat  diameter)  from  639  down  to  57. 

Two  of  the  nozzle  flows  presented  in  Ref  10  were  chosen  as  test  cases  for 
the  present  work.  In  both  examples  the  'propellant'  is  nitrogen  at  an  initial 
temperature  of  300  K.  Also  in  both  examples  the  semi-angle  of  the  convergent 
part  of  the  nozzle  is  30°,  the  semi-angle  of  the  nozzle  expansion  cone  is  20°  and 
the  longitudinal  radius  of  curvature  of  the  throat  is  one  half  the  transverse 
radius.  For  the  purposes  of  this  computation  the  radius  of  the  inlet  was  taken 
to  be  three  times  the  throat  radius. 

5. 1 Example  1 

The  throat  Reynolds  number  (based  on  diameter)  = 639  and  the  throat 

3 2 -2 

diameter  * 2.50  mm,  the  inlet  pressure  *=  1.981  x 10  N/m  (1.955  x io  bar)  and 
the  mass  flow  rate  = 20.1  mg/s.  The  computed  flow-field  is  shown  in  Fig  3.  The 
lower  half  of  the  nozzle  shows  the  streamtubes  and  orthogonal  (constant  pressure) 
surfaces.  The  upper  half  shows  the  displacement  thickness  and  constant  Mach 
surfaces  at  M * 1,  2,  3 and  4.  The  computation  was  performed  with  20  streamtubes 
although  only  10  are  shown  for  clarity.  175 
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It  will  be  noted  that  the  displacement  thickness  is  compressed  consider- 
ably by  the  convergent  part  of  the  nozzle,  but  grows  out  into  the  stream  in  the 
expansion  cone.  Although  about  95%  of  the  mass  flow  becomes  sonic  only  about 
70%  is  accelerated  to  Mach  3. 

It  is  unsafe  to  make  conclusions  of  an  analytic  nature  from  a numerical 
analysis,  especially  in  the  region  of  singularities,  but  the  behaviour  of  the 
sonic  surface  at  the  nozzle  throat  may  be  compared  with  the  sonic  surface  in 
two-dimensional  inviscid  flow.  In  inviscid  flow  the  sonic  surface  is  parabolic 
in  form  near  the  symmetry  axis,  intersects  the  axis  downstream  of  the  geometric 
throat  and  is  concave  towards  the  flow  (see  for  example  Ref  19).  Here  the  sonic 
surface  is  slightly  convex  towards  the  flow. 

The  axial  static  temperature  is  shown  in  Fig  4.  Near  the  throat  the  agree- 
ment with  experimental  results  is  excellent  and  even  near  the  exit  plane  the 
computed  curve  is  only  just  outside  the  error  bars.  Fig  5 shows  the  results  for 
the  axial  static  pressure.  The  computed  curve  agrees  with  the  experimental 
curve  to  within  20%. 

5.2  Example  2 

The  throat  Reynolds  number  - 307  and  the  throat  diameter  = 5.1  mm,  the 

2 2 -3 

inlet  pressure  *■  4.657  x 10  N/m  (4.596  x 10  bar)  and  the  mass  flow  rate 
18.97  mg/s. 

The  computed  flow-field  is  shown  in  Fig  6.  The  displacement  thickness 
curves  more  rapidly  into  the  flow  in  the  expansion  cone  and  only  about  45%  of 
the  mass  flow  is  accelerated  to  M * 3. 

The  axial  static  temperature  is  plotted  in  Fig  7 and  here  again  the  agree- 
ment with  experimental  results  is  good.  Note  the  oscillation  in  the  curve  at 
X/Rt  * 17  due  to  an  incipient  instability.  Fig  8 shows  the  axial  static 
pressure  and  here  the  agreement  is  excellent  except  for  the  oscillation  near 
X/Rt  ■ 17.  The  cause  of  this  oscillation  has  not  yet  been  identified.  A further 
example  at  a lower  Reynolds  number  (91)  was  tried,  but  increased  instability 
prevented  a reliable  computation. 

In  the  present  case  Rothe  presents  results  for  the  radial  profiles  of 
temperature,  density  and  pressure  at  X/Rt  - 13.7.  Fig  9 shows  the  radial  static 
temperature.  The  computed  results  are  within  or  on  the  edge  of  the  error  bars 
over  the  whole  range.  In  Fig  10  the  radial  density  profiles  are  shown,  here  the 
agreement  is  not  so  good  but  still  reasonable.  Fig  1 I shows  the  radial  pressure 
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profile  and  here  the  agreement  is  poor.  At  the  nozzle  wall  the  computed 
pressure  loss  is  about  121,  but  the  measured  loss  is  about  29. 5%.  Fig  6 shows 
that  the  streamtubes  are  only  very  slightly  curved  in  this  region  so  that  the 
extra  pressure  loss  is  unlikely  to  be  caused  by  a centrepetal  force  component, 
and  in  any  event  the  streamtubes  are  curved  towards  the  symmetry  axis  which  would 
tend  to  increase  the  pressure  along  the  radius.  It  is  tempting  to  speculate  that 
the  pressure  effect  is  caused  by  the  influence  of  the  low  pressure  in  the  vacuum 
chamber  in  which  the  nozzle  is  tested,  being  propagated  upstream  through  the 
subsonic  part  of  the  flow.  The  procedure  described  in  this  Report,  being  based 
on  a parabolic  approximation,  is  unable  to  take  this  kind  of  effect  into  account. 

6 CONCLUSIONS 

The  object  of  the  work  presented  in  this  Report  was  to  develop  a numerical 
procedure  suitable  for  computing  the  entire  flow-field  in  a nozzle  at  low 
Reynolds  numbers.  As  the  results  in  examples  I and  2 show  this  has  been 
reasonably  successful.  No  approximations  are  required  other  than  those  inherent 
in  the  boundary  layer  equations.  Similarity  is  not  assumed.  A viscous  boundary 
layer- inviscid  cone  flow  is  not  assumed.  Furthermore  the  method  is  of  the 
direct  type,  the  flow  can  be  computed  from  a given  nozzle  geometry. 

Despite  the  low  density  of  the  flow  in  the  examples,  non-slip  boundary 
conditions  were  employed  and  gave  reasonable  results. 

A criticism  sometimes  levelled  at  full  flow-field  computations  is  that 
they  are  slow  and  require  large  fast  computers.  This  may  have  been  true  in  the 
past,  but  computers  are  becoming  more  powerful,  so  much  so  that  procedures  of  the 
type  described  here  are  now  practicable  propositions. 
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NOMENCLATURE 

A arbitrary  coefficient  in  equation  (30),  defined  in  identities  31  and  34, 
or  area  in  equation  (60) 
a radius  of  nozzle  wall 

Sq  radius  of  nozzle  inlet 

a(  radius  of  nozzle  throat 

' 

B arbitrary  coefficient  in  equations  (28)  and  (30),  defined  in  identities  32 
and  35 

b defined  in  equation  (42) 

C arbitrary  coefficient  in  equation  (30),  defined  in  identities  33  and  36 
Cp  specific  heat  at  constant  pressure 

j!  cp  cp/cpo 

c matrix  coefficient  defined  in  equations  (38),  (39)  and  (40) 
h static  enthalpy  (Ah  - CpAT) 

h*  h/hQ 

k thermal  conductivity 

k*  k/kQ 

M Mach  number 

m molecular  weight 

m mass  flow  rate 


N number  of  streamtubes 

n et'afTr  nrpecurA 


20 


v 

X 

x 

X* 

y 

y* 

Y 

A 

e 

» 


x/«t  , 


NOMENCLATURE  (concluded) 

velocity  in  transverse  direction 

x measured  along  symmetry  axis  from  nozzle  throat 

coordinate  in  longitudinal  direction 
(x/a0)/proReo 

coordinate  in  transverse  direction 

y/*o 

ratio  of  specific  heats 
finite  difference 
streamtube  inclination 
viscosity 
u/u. 


0 

p density 

P*  p/pQ 

♦ arbitrary  function 

<|i  stream  function 

Subscripts 

m denotes  arbitrary  point  in  x direction 

mean  denotes  mean  value 

n denotes  arbitrary  point  in  direction 

0 denotes  initial  conditions 

v denotes  value  at  nozzle  vail 

Superscripts 

* denotes  non-dimensional  variable 
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Figs  1&2 


Fig  1 The  coordinate  system 


Fig  2 General  point  in  the  mesh  scheme  and  domain  of  integration  x*— plane 
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17.  Abstract 


A numerical  procedure  for  computing  the  laminar  flow-field  in  nozzles  at 
throat  Reynolds  numbers  of  300-3000  is  described.  Such  nozzles  are  found  in  space- 
craft position  control  thrusters,  chemical  lasers,  and  low  density  hypersonic  wind- 
tunnels.  A parabolic  approximation  of  the  Navier-Stokes  equations  (the  /^boundary 
layer  equations*;  is  transformed  to  Von  Mises  form  and  solved  numerically  by  a 
central  difference  method.  The  entire  subsonic  and  supersonic  flow-field  is 
computed.  No  assumptions  or  approximations,  other  than  those  inherent  in  the  flow 
equations,  are  involved.  The  method  is  of  the  direct  type  and  the  flow-field  for 
any  given  nozzle  geometry  may,  in  principle,  be  computed.  Examples  are  given 
comparing  computed  results  with  published  experimental  data. 
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